## ---------------------------------------------------------------------- #
## Figure 8
## Christopher J. Fariss
## "Yes, Human Rights Practices Are Improving Over Time"
## American Political Science Review
## ---------------------------------------------------------------------- #
rm(list = ls())

#7 binary event variables: genocide rummel massive_repression executions negative_sanctions mass_killing killing_prsent

#label these EVENTS


#Model1: EVENTS
#Model2: EVENTS + HRW
#Model3: EVENTS + HRW + ITT
#Model4: EVENTS + HRW + ITT + POLPRIS
#Model5: EVENTS + HRW + ITT + POLPRIS + DISAP
#Model6: EVENTS + HRW + ITT + POLPRIS + DISAP + Amnesty
#Model7: EVENTS + HRW + ITT + POLPRIS + DISAP + Amnesty + KILL
#Model8: EVENTS + HRW + ITT + POLPRIS + DISAP + Amnesty + KILL + hathaway
#Model9: EVENTS + HRW + ITT + POLPRIS + DISAP + Amnesty + KILL + hathaway + TORT
#Model10: EVENTS + HRW + ITT + POLPRIS + DISAP + Amnesty + KILL + hathaway + TORT + State


year <- 1946:2015



NAMES <- c("Event-based items only (7 items)",
"+ PTS: Human Rights Watch (8 items)",
"+ ITT: Torture (9 items)",
"+ CIRI: Political Imprisonment (10 items)",
"+ CIRI: Disappearance (11 items)",
"+ PTS: Amnesty (12 items)",
"+ CIRI: Extrajudicial Killing (13 items)",
"+ Hathaway: Torture (14 items)",
"+ CIRI: Torture (15 items)" ,
"+ PTS: State (16 items)"
)


for(i in 1:9){
    
    MODEL <- i
    
    
    if(MODEL==1) load("progressive_summary_events.RData")
    if(MODEL==2) load("progressive_summary_hrw.RData")
    if(MODEL==3) load("progressive_summary_itt.RData")
    if(MODEL==4) load("progressive_summary_polpris.RData")
    if(MODEL==5) load("progressive_summary_disap.RData")
    if(MODEL==6) load("progressive_summary_amnesty.RData")
    if(MODEL==7) load("progressive_summary_kill.RData")
    if(MODEL==8) load("progressive_summary_hathaway.RData")
    if(MODEL==9) load("progressive_summary_tort.RData")
    
    
    if(MODEL==1) test <- read.csv("progressive_summary_events.csv")
    if(MODEL==2) test <- read.csv("progressive_summary_hrw.csv")
    if(MODEL==3) test <- read.csv("progressive_summary_itt.csv")
    if(MODEL==4) test <- read.csv("progressive_summary_polpris.csv")
    if(MODEL==5) test <- read.csv("progressive_summary_disap.csv")
    if(MODEL==6) test <- read.csv("progressive_summary_amnesty.csv")
    if(MODEL==7) test <- read.csv("progressive_summary_kill.csv")
    if(MODEL==8) test <- read.csv("progressive_summary_hathaway.csv")
    if(MODEL==9) test <- read.csv("progressive_summary_tort.csv")
    
    
    
    tmp <- subset(test, grepl("theta", test$X), select=c(mean, sd))
    
    
    data  <- data[order(data$id),]
    
    data$latentmean <- tmp$mean
    data$latentsd <- tmp$sd
    
    year <- YEAR <- 1946:2015
    
    temp1 <- temp2 <- NA
    for(i in 1:length(year)){
        m <- data$latentmean[data$YEAR==year[i]]
        v <- data$latentsd[data$YEAR==year[i]]
        temp <- apply(matrix(rnorm(matrix(1,ncol=1000, nrow=length(m)), mean=m, sd=v),ncol=1000, nrow=length(m)), 2, mean)
        temp1[i] <- mean(temp)
        temp2[i] <- sd(temp)
    }
    latentmean <- temp1
    latentsd <- temp2
    
    data <- as.data.frame(cbind(YEAR, latentmean, latentsd))
    
    COLOR <- c(2,3)
    COLOR <- c(1,1)
    
    temp <- matrix(-999, nrow=5000, ncol=length(year))
    
    # take 5000 draws from a normal distribution for each year of the data, which fills 5000 rows in each i column of the data
    for(i in year-1945){
        temp[,i]  <- rnorm(5000, mean=data$latentmean[data$YEAR==year[i]], sd=data$latentsd[data$YEAR==year[i]])
    }
    
    #par(mar=c(2,3,0.2,2), cex=1.2, font=2, font.lab=2)
    par(mar=c(2,5,2,2), cex=1.2, font=1, font.lab=1)
    
    # calculate the mean value of each column in the matrix; i.e., the mean of each column of 5000 draws
    point.estimates <- apply(temp, 2, mean)
    
    # plot the country-year point.estimates
    plot(point.estimates, ylim=c(-0.6,1.20), yaxt="n", xaxt="n", pch=19, col=grey(.25), xlab="", ylab="", cex=.8, type="n", xlim=c(1,70))
    axis(1, at=(1:70), label=rep("", length(1:70)))
    axis(1, at=(seq(5,65,10)), label=seq(1950,2010,10), lwd=1.5)
    axis(2, at=c(-0.75,-0.50,-.25,0,0.25,0.50,0.75, 1.0), las=2)
    #axis(2, at=c(-1,-0.5,0,0.5,1.0), las=2)
    mtext(side=2, "Latent Physical Integrity", line=3.5, cex=1.7)
    mtext(side=4, "More Abuse", font=1, line=0.3, cex=1.7, at=-0.15)
    mtext(side=4, "More Respect", font=1, line=0.3, cex=1.7, at=0.65)
    arrows(75.25, -.50, 75.25, -.65, xpd = TRUE, lwd=3, length=.15)
    arrows(75.25, 1.05, 75.25, 1.20, xpd = TRUE, lwd=3, length=.15)
    
    abline(h=-.25,lty=2, col=grey(.75), lwd=2)
    abline(h=0,lty=2, col=grey(.75), lwd=2)
    abline(h=.25,lty=2, col=grey(.75), lwd=2)
    abline(h=.5,lty=2, col=grey(.75), lwd=2)
    abline(h=.75,lty=2, col=grey(.75), lwd=2)
    abline(h=1.00,lty=2, col=grey(.75), lwd=2)
    
    # add 95% credible intervales to the plot
    for(i in 1:length(year)){
        lines(c(i,i), c(apply(temp, 2, quantile, .025)[i], apply(temp, 2, quantile, .975)[i]), lwd=2, col=grey(.6))
    }
    points(apply(temp, 2, mean), col=grey(.6), pch=21, bg=grey(.95), cex=.8)
    
    
    
    
    mtext(side=3, NAMES[MODEL], line=0.3, cex=1.7)
}




data <- read.csv("M_1_Full_Data.csv")
data$latentmean <- data$theta_mean
data$latentsd <-data$theta_sd

year <- YEAR <- 1946:2015

temp1 <- temp2 <- NA
for(i in 1:length(year)){
    m <- data$latentmean[data$YEAR==year[i]]
    v <- data$latentsd[data$YEAR==year[i]]
    temp <- apply(matrix(rnorm(matrix(1,ncol=1000, nrow=length(m)), mean=m, sd=v),ncol=1000, nrow=length(m)), 2, mean)
    temp1[i] <- mean(temp)
    temp2[i] <- sd(temp)
}
latentmean <- temp1
latentsd <- temp2

data <- as.data.frame(cbind(YEAR, latentmean, latentsd))

COLOR <- c(2,3)
COLOR <- c(1,1)

temp <- matrix(-999, nrow=5000, ncol=length(year))

# take 5000 draws from a normal distribution for each year of the data, which fills 5000 rows in each i column of the data
for(i in year-1945){
    temp[,i]  <- rnorm(5000, mean=data$latentmean[data$YEAR==year[i]], sd=data$latentsd[data$YEAR==year[i]])
}


# calculate the mean value of each column in the matrix; i.e., the mean of each column of 5000 draws
point.estimates <- apply(temp, 2, mean)

# plot the country-year point.estimates
plot(point.estimates, ylim=c(-0.6,1.20), yaxt="n", xaxt="n", pch=19, col=grey(.25), xlab="", ylab="", cex=.8, type="n", xlim=c(1,70))
axis(1, at=(1:70), label=rep("", length(1:70)))
axis(1, at=(seq(5,65,10)), label=seq(1950,2010,10), lwd=1.5)
axis(2, at=c(-0.75,-0.50,-.25,0,0.25,0.50,0.75, 1.0), las=2)
#axis(2, at=c(-1,-0.5,0,0.5,1.0), las=2)
mtext(side=2, "Latent Physical Integrity", line=3.5, cex=1.7)
mtext(side=4, "More Abuse", font=1, line=0.3, cex=1.7, at=-0.15)
mtext(side=4, "More Respect", font=1, line=0.3, cex=1.7, at=0.65)
arrows(75.25, -.50, 75.25, -.65, xpd = TRUE, lwd=3, length=.15)
arrows(75.25, 1.05, 75.25, 1.20, xpd = TRUE, lwd=3, length=.15)

abline(h=-.25,lty=2, col=grey(.75), lwd=2)
abline(h=0,lty=2, col=grey(.75), lwd=2)
abline(h=.25,lty=2, col=grey(.75), lwd=2)
abline(h=.5,lty=2, col=grey(.75), lwd=2)
abline(h=.75,lty=2, col=grey(.75), lwd=2)
abline(h=1.00,lty=2, col=grey(.75), lwd=2)

# add 95% credible intervales to the plot
for(i in 1:length(year)){
    lines(c(i,i), c(apply(temp, 2, quantile, .025)[i], apply(temp, 2, quantile, .975)[i]), lwd=2, col=grey(.6))
}
points(apply(temp, 2, mean), col=grey(.6), pch=21, bg=grey(.95), cex=.8)

mtext(side=3, NAMES[10], line=0.3, cex=1.7)


